;; Tests of the Non-robust computational geometry functions
;; Copyright (C) 2022  Rocks Mazama
;; 
;; This program is free software: you can redistribute it and/or modify
;; it under the terms of the GNU General Public License as published by
;; the Free Software Foundation, either version 3 of the License, or
;; (at your option) any later version.
;; 
;; This program is distributed in the hope that it will be useful,
;; but WITHOUT ANY WARRANTY; without even the implied warranty of
;; MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
;; GNU General Public License for more details.
;; 
;; You should have received a copy of the GNU General Public License
;; along with this program.  If not, see <https://www.gnu.org/licenses/>.

(ns mazama.compgeom.core-test
  (:require
   [clojure.test :refer [deftest is testing]]
   [mazama.compgeom.core :refer :all]))

;;;; NOTICE
;;;;
;;;; It is assumed that
;;;; the x-axis increases to the right and
;;;; the y-axis increases DOWNWARD

(def EPSILON 1/1000000)
 
;;; Reference: https://floating-point-gui.de/errors/comparison/
(defn nearly==
  ""
  ([a b] (nearly== a b EPSILON))
  ([a b epsilon]
   (or (== a b)
       (if (and (rational? a) (rational? b))
         false
         (let [absA (Math/abs a)
               absB (Math/abs b)
               diff (Math/abs (- a b))]
           (if (or (== a 0) (== b 0) (< (+ absA absB) Float/MIN_NORMAL))
             (< diff (* epsilon Float/MIN_NORMAL))
             (< (/ diff (min (+ absA absB) Float/MAX_VALUE)) epsilon)))))))

(defn nearly<=
  ""
  ([a b] (nearly<= a b EPSILON))
  ([a b epsilon] (or (<= a b) (nearly== a b epsilon))))

(defn nearly>=
  ""
  ([a b] (nearly>= a b EPSILON))
  ([a b epsilon] (or (>= a b) (nearly== a b epsilon))))

(defn perturb [dx dy [vx vy]]
  [(+ vx dx) (+ vy dy)])

(defn scale [t [vx vy]]
  [(* t vx) (* t vy)])

(deftest foot-test
  (testing "colinear"
    (let [A [-1 2]
          V [-3 -2]
          P [2 4]
          Q A
          R [-5/2 1]
          S [-4 0]
          T [-7 -2]]
      (is (= -1 (foot P A V)))
      (is (= 0 (foot Q A V)))
      (is (= 1/2 (foot R A V)))
      (is (= 1 (foot S A V)))
      (is (= 2 (foot T A V)))))
  (testing "horizontal"
    (let [A [3 -1]
          V [4 0]
          P [2 10]
          Q [4 -2]
          R [9 -4]
          S [11 -1]]
      (is (= -1/4 (foot P A V)))
      (is (= 1/4 (foot Q A V)))
      (is (= 3/2 (foot R A V)))
      (is (= 2 (foot S A V)))))
  (testing "vertical"
    (let [A [2 1]
          V [0 -3]
          P [0 3]
          Q [3 0]
          R [5 -3]
          S [1 7]]
      (is (= -2/3 (foot P A V)))
      (is (= 1/3 (foot Q A V)))
      (is (= 4/3 (foot R A V)))
      (is (= -2 (foot S A V))))))

(deftest foot-distance-test
  (testing "horizontal"
    (is (==  1 (foot-distance [-1 -1] [-1/2 0] [1 0])))
    (is (==  1 (foot-distance [ 0 -1] [-1/2 0] [1 0])))
    (is (==  1 (foot-distance [ 1 -1] [-1/2 0] [1 0])))
    (is (== -1 (foot-distance [-1  1] [-1/2 0] [1 0])))
    (is (== -1 (foot-distance [ 0  1] [-1/2 0] [1 0])))
    (is (== -1 (foot-distance [ 1  1] [-1/2 0] [1 0]))))
  (testing "vertical"
    (is (== -1 (foot-distance [ 1 -1] [0 1/2] [0 -1])))
    (is (== -1 (foot-distance [ 1  0] [0 1/2] [0 -1])))
    (is (== -1 (foot-distance [ 1  1] [0 1/2] [0 -1])))
    (is (==  1 (foot-distance [-1 -1] [0 1/2] [0 -1])))
    (is (==  1 (foot-distance [-1  0] [0 1/2] [0 -1])))
    (is (==  1 (foot-distance [-1  1] [0 1/2] [0 -1]))))
  (testing "diagonal"
    (is (nearly== (/ (Math/sqrt 2) 2) (foot-distance [1 0] [0 0] [1 1]))))
  (testing "Example 1: https://brilliant.org/wiki/dot-product-distance-between-point-and-a-line/"
    (is (nearly== (/ 3 (* 2 (Math/sqrt 5))) (foot-distance [-3 2] [0 5/4] [1 -1/2]))))
  (testing "Example 2: https://brilliant.org/wiki/dot-product-distance-between-point-and-a-line/"
    (is (== 6/5 (foot-distance [0 0] [0 3/2] [1 -3/4])))))

(deftest parametric-line-test
  (testing "input parameter - return the point"
    (let [a [3 -1]
          b [7 9]
          eq (parametric-line a b)]
      (is (= a (eq 0)))
      (is (= b (eq 1)))
      (is (= [5 4] (eq 1/2)))))
  (testing "input point - return the parameter of the foot"
    (let [a [0 0]
          b [1 3]
          eq (parametric-line a b)]
      (is (= 1/2 (eq 7/2 1/2))))))

(deftest point-orient2d-test
  (testing "ccw"
    (is (= 1 (point-orient2d [0 0] [0 1] [1 1])))
    (is (= 1 (point-orient2d [1 0] [11 11] [22 22] )))
    (is (= 1 (point-orient2d [5 5] [3 3] [99 100]))))
  (testing "cw"
    (is (= -1 (point-orient2d [-3 -2] [1 1] [-1 2])))
    (is (= -1 (point-orient2d [123456 123456] [0 0] [123457 123455])))
    (is (= -1 (point-orient2d [0 0] [12345 12344] [123457 123457]))))
  (testing "colinear"
    (is (= 0 (point-orient2d [-4 -4] [-4 -4] [-4 -4])))
    (is (= 0 (point-orient2d [1 1] [1 1] [10 100])))
    (is (= 0 (point-orient2d [2 2] [-1 10] [2 2])))
    (is (= 0 (point-orient2d [3 7] [2 2] [2 2])))
    (is (= 0 (point-orient2d [999999 -999999] [1 -1] [-98765 98765])))
    (is (= 0 (point-orient2d [-11/10 -11/10] [99/10 99/10] [33/10 33/10])))
    (is (= 0 (point-orient2d [333/10 333/10] [999/10 999/10] [-111/10 -111/10])))
    ))

(deftest vector-orient2d-test
  (testing "ccw"
    (is (= 1 (vector-orient2d [1 1] [1 0])))
    (let [v [12345 67890]
          w (perturb 0 (- EPSILON) (scale 3 v))
          result (vector-orient2d v w)]
      (is (= 1 result)))
    (let [v [12345 67890]
          w (perturb EPSILON 0 (scale 1/3 v))
          result (vector-orient2d v w)]
      (is (= 1 result)))
    (let [v [9876 54321]
          w (perturb 0 (- EPSILON) (scale -4 v))
          result (vector-orient2d v w)]
      (is (= 1 result)))
    (let [v [9876 54321]
          w (perturb EPSILON 0 (scale -2/5 v))
          result (vector-orient2d v w)]
      (is (= 1 result))))
  (testing "cw"
    (is (= -1 (vector-orient2d [1 0] [1 1])))
    (let [v [12345 67890]
          w (perturb 0 EPSILON (scale 3 v))
          result (vector-orient2d v w)]
      (is (= -1 result)))
    (let [v [12345 67890]
          w (perturb (- EPSILON) 0 (scale 3/10 v))
          result (vector-orient2d v w)]
      (is (= -1 result)))
    (let [v [9876 54321]
          w (perturb 0 EPSILON (scale -4 v))
          result (vector-orient2d v w)]
      (is (= -1 result)))
    (let [v [9876 54321]
          w (perturb (- EPSILON) 0 (scale -2/5 v))
          result (vector-orient2d v w)]
      (is (= -1 result))))
  (testing "colinear"
    (is (= 0 (vector-orient2d [0 0] [1 1])))
    (is (= 0 (vector-orient2d [-1 -1] [0 0])))
    (let [v [2 4]
          w (scale 2/3 v)
          result (vector-orient2d v w)]
      (is (= 0 result)))
    (let [v [-2 3]
          w (scale 22/7 v)
          result (vector-orient2d v w)]
      (is (= 0 result)))
    (let [v [-7 -6]
          w (scale 13/11 v)
          result (vector-orient2d v w)]
      (is (= 0 result)))
    (let [v [2/3 -22/7]
          w (scale -37/3 v)
          result (vector-orient2d v w)]
      (is (= 0 result)))
    ))

(deftest vector-angle-test
  (testing "degenerate"
    (is (== 0.0 (vector-angle [0 0] [0 0])))
    (is (== 0.0 (vector-angle [0 0] [2 3])))
    (is (== 0.0 (vector-angle [-4 -5] [0 0]))))
  (is (== (* Math/PI -3/4) (vector-angle [1 0] [-1 1])))
  (is (== (/ Math/PI 2) (vector-angle [0 1] [1 0])))  
  (is (== Math/PI (vector-angle [1.1 1.1] [-2.2 -2.2]))))

(deftest interval-intersect?-test
  (testing "false"
    (is (= false (interval-intersect? [0 1] [2 3])))
    (is (= false (interval-intersect? [2 3] [0 1])))
    (is (= false (interval-intersect? [2 2] [3 3])))
    (is (= false (interval-intersect? [2 3] [31/10 4]))))
  (testing "true"
    (is (= true (interval-intersect? [2 2] [2 2])))
    (is (= true (interval-intersect? [5 8] [4 5])))
    (is (= true (interval-intersect? [5 5] [5 8])))
    (is (= true (interval-intersect? [5 8] [4 6])))
    (is (= true (interval-intersect? [6 7] [5 8])))
    (is (= true (interval-intersect? [5 8] [8 8])))
    (is (= true (interval-intersect? [7 9] [5 8])))))

(deftest incircle-test
  (testing "ccw-inside"
    (is (pos? (incircle [-1 0] [0 1] [1 0] [0 0]))))
  (testing "cw-inside"
    (is (neg? (incircle [3 1] [1 3] [-1 1] [1 1]))))
  (testing "ccw-cocircular"
    (is (zero? (incircle [-1 -1] [0 0] [1 -1] [0 -2]))))
  (testing "cw-cocircular"
    (is (zero? (incircle [0 1] [-1 0] [0 -1] [1 0]))))
  (testing "ccw-outside"
    (is (neg? (incircle [0 -1] [-1 0] [0 1] [-241/239 0]))))
  (testing "cw-outside"
    (is (pos? (incircle [1 0] [0 1] [-1 0] [0 -2])))))

(deftest general-line-test
  (let [coll [[[0 0] [1 1]]
              [[1 2] [-2 1]]
              [[-1/2 1/3] [2/5 -5/2]]]]
    (doseq [[[x y :as p] [vx vy :as v]] coll]
      (let [[a b c] (general-line p v)]
        (doseq [t [0 1 -1 3/10 10/3 -10/3 1/2]]
          (is (zero? (+ (* a (+ x (* t vx))) (* b (+ y (* t vy))) c))))))))

(deftest line-intersection-test
  (is (= [-6/5 7/5] (line-intersection [3 4 -2] [2 1 1])))
  (is (= [-3 2] (line-intersection [-1 3 -9] [2 1 4])))
  (is (= [1 -2] (line-intersection [2 4 6] [2 3 4])))
  (is (= [-7 3] (line-intersection [2 4 2] [2 3 5])))
  (is (= [10 27] (line-intersection [-3 1 3] [-23/10 1 -4])))
  (is (= [-2/5 11/5] (line-intersection [2 -1 3] [1 2 -4])))
  (is (= [-20/17 46/17] (line-intersection [4 1 2] [-1/4 1 -3]))))
